Shear flow of non-Brownian suspensions close to jamming 
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The dynamical mechanisms controlling the rheology of dense suspensions close to jamming are 
investigated numerically, using simplified models for the relevant dissipative forces. We show that 
the velocity fluctuations control the dissipation rate and therefore the effective viscosity of the 
suspension. These fluctuations are similar in quasi-static simulations and for finite strain rate 
calculations with various damping schemes. We conclude that the statistical properties of grain 
trajectories - in particular the critical exponent of velocity fluctuations with respect to volume 
fraction cj> - only weakly depend on the dissipation mechanism. Rather they are determined by 
steric effects, which are the main driving forces in the quasistatic simulations. The critical exponent 
of the suspension viscosity with respect to <j> can then be deduced, and is consistent with experimental 
data. 

PACS numbers: 66.20. Cy,83.80.Hj 



Athermal disordered systems such as foams [l|, emul- 
sions [2 L non-Brownian suspensions Q or granular ma- 
terials 4J exhibit a critical phase transition between a 
liquid-like and a solid-like mechanical behaviour, when 
the particle volume fraction (f> crosses the jamming point 
4> c . For (f> > <p c , these amorphous systems can resist 
shear. The elastic shear modulus vanishes at <f> c with a 
critical exponent different from the mean field one 0] . 
Above a yield stress ay, vanishing at (f> Cl they present a 
non-Newtonian rheology, for which several different inter- 
pretations have been proposed, based on (i) an analogy 
with the glassy dynamics of a system presenting scale- free 
energy distributions Q , (ii) interacting plastic events [|| , 
(iii) the critical scaling laws of the shear modulus and of 
the coordination number 0] . Together with conventional 
molecular dynamics simulations (MD), quasistatic meth- 
ods (QS) have been applied to study the plastic flow of 
athermal amorphous solids at the yield-stress ay 
It is generally assumed that QS accurately describe the 
dynamics of the true system in the limit of asymptotically 
small shear rate 7. However, the existence of a proper 
quasistatic limit remains controversial, and there is grow- 
ing evidence that quasistatic flows actually correspond to 
a finite-size dominated regime, with a correlation length 
that saturates at the system size [1, [3] ■ 

Symmetrically, for <f> < <p c , amorphous materials can 
flow under an infinitesimal shear stress a and present a 
viscosity 77 diverging at <f> c like r\ cx (<p c — <fi)~~ a ■ Scal- 
ing laws are expected to be different in thermal (glassy) 
systems and in athermal sytems [l6| . In the case of a sus- 
pension of non-Brownian particles, the best fit of recent 
experimental results give a critical exponent of a = 2.4 
for volume-controlled experiments [HI and of a = 1.9 for 
pressure-controlled experiments 0]. The explanation of 
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FIG. 1. (Color online) Viscosity 77 normalized by the friction 
coefficient £ as a function of the volume fraction <f>, measured 
from molecular dynamics (MD) and quasi static simulations 
(QS), for a dissipation induced by viscous drag forces of Eq. 
[TJ (labelled v — 1) or by the lubrication like mechanism of Eq. 
(labelled LF). In MD, the viscosity is measured in the low 
shear rate regime, for 7 = 10 -6 and £ = 10 _1 . In this regime, 
n/( does not depend on the precise value of these parameters, 
as shown in Fig. [2k. The quantitative agreement between 
MD(z/ = 1) and MD(LF) is coincidental. Inset: compilation 
of experimental data available in the literature at imposed 
pressure P (with 4> c = 0.587) and imposed volume fraction 
<j> (with <j> c = 0.615) for the ratio of suspension to solvent 
viscosity. 
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the critical exponent as well as the underlying mecha- 
nisms of the flow arrest have remained open and con- 
troversial questions up to now. Among the proposed 
mechanisms are hydrodynamic dissipation in the lubri- 
cated films separating particles, or friction-induced nor- 
mal stresses [17] . A completely different interpretation 
relates the divergence of the viscosity to a singular mode 
of the network of contacts close to the isostatic point [3] ■ 

Here, we present simulation results for the viscous flow 
of a simplified model system in the vicinity of the close- 
packed state at <f> c . We identify a dynamical contribution 
to the divergence of the viscosity, which has its origin 
in the singularity of velocity fluctuations. By compar- 
ing different computational model systems we show that 
some of the statistical properties of these velocity fluc- 
tuations are surprisingly model independent. We show 
how the rheological properties, in particular the form of 
the flow curve and the divergence of the viscosity can be 
obtained from one set of trajectories that is based on a 
quasistatic simulation method. 

Simulations - We consider a two-dimensional system 
constituted by N soft spherical particles of mass m, N/2 
of diameter d and N/2 of diameter l.Ad. The particle vol- 
ume fraction is defined as </> = Yli=i 7rr ?/^ 2 j where L is 
the size of the simulation box. Periodic (Lees-Edwards) 
boundary conditions are used in both directions. Two 
particles i,j interact when their distance r is smaller than 
the sum of their radii r, + rj , with a repulsive potential 
E(r) = e(l — r/(rj + rj)) 2 . All observables below are 
given in units of d, m and e. We compare the divergence 
of the viscosity for tp < (f> c = 0.843 ID, lj| using two dif- 
ferent dynamics: non-equilibrium dissipative molecular 
dynamics (MD) and quasistatic simulations (QS). 

In the MD simulations, the system is sheared at a shear 
rate 7. Newton's equations of motion mri = F? + i<7 ISC 
are integrated with elastic contact forces F cl = —VE and 
a viscous drag force 



(1) 



proportional to the isth power of the velocity difference 
5vi = Vi — ?fl ow between the particle velocity Vi and the 
flow velocity i?flnw(^.) = Zxiy-, whose fluctuations are ne- 
glected 23-23- The flow can be viewed as being set 
up by a non-Newtonian fluid characterized by a friction 
coefficient (. In the special case v = 1, the fluid is 
Newtonian and C is proportional to the bare fluid vis- 
cosity rjf (in Stokes approximation, £ = i-Krjfd). Ther- 
mal and lubrication forces are ignored. Unlike in gran- 
ular systems, the particle-particle collisions are elastic 
and the only dissipation is due to viscous losses associ- 
ated with the fluctuations of the particle velocity field. 
The shear stress a is calculated from the particle posi- 
tions fi = (xi,yi) and the forces Fi = (Fi Xl Fi y ) acting 
on them as a = L^ 2 ^2^ =1 XiFi y . The dominant contri- 
bution comes from the elastic forces that result from par- 



ticle overlaps. The resulting relation between the shear 
stress a and the shear rate 7 is shown in Fig. [2^. For 
small strain rates, both inertia and deformation of the 
particles are negligible, and the stress grows with strain 
rate as a = iff" ', characteristic for a power-law fluid. The 
"effective viscosity" r](4>) is measured in this regime and 
is a function of the volume-fraction, as shown in Fig. [1] 
for the Newtonian case (v = 1). At larger strainrates and 
in weakly damped systems (£ = 0.001 or v > 1) one ob- 
serves a shear thickening regime, which can be ascribed 
to inertia. Conversly, for stronger damping (f = 0.1 and 
v < 1) and, in particular for volume fractions close to <j) c 
[23] one observes a shear thinning regime, when particle 
deformation starts to be relevant. 

Quasistatic simulations consist of successively applying 
small steps of shear and minimizing the total potential 
energy. By construction, they generate particle trajecto- 
ries at 7 0. An elementary strain step of 70 = 5 ■ 10~ 5 
is used. After each change in boundary conditions the 
particles are moved affincly to define the starting config- 
uration for the minimization, which is performed using 
conjugate gradient techniques [24j |. The minimization is 
stopped when the nearest energy minimum is found. As 
no static, force-balanced state exists below the jamming 
point ((f) < (f> c ), the inter-particle forces at the minimum 
are strictly zero; i.e. the particles can always arrange in 
such a way as to avoid mutual overlaps. Thus, each min- 
imized configuration corresponds to a true hard-sphere 
state and the resulting particle trajectories can be viewed 
as a sequence of snapshots of a flowing hard-sphere sys- 
tem at zero temperature. Particle motion in such a sys- 
tem is driven by steric exclusion and the lack of free vol- 
ume. In particular, particles have to move over larger 
distances when the jamming point is approached, to find 
a new overlap-free state compatible with the imposed 
shear j2j|. 

Without particle overlaps all contact forces and there- 
fore the shear stress are strictly zero in the QS simula- 
tion. Still, an effective shear stress and viscosity can be 
obtained from the power T per unit surface that would 
be dissipated along the QS trajectories, if the dissipa- 
tion mechanism of Eq. ([T]) was present. T is equal to the 
power injected per unit volume in the system, cry, and 
can be expressed as: 

From this expression, we deduce the viscosity: 



T N (5v 1+u ) N f ,, 

" = 7^ = -<Z» V 7 ^ = ~ C T 2 J Al+ ^( A ) dA ■ 

(2) 

where -P(A) is the probability distribution function of the 
particle velocity rescaled by the shear rate: A = Sv/'j. 
As particle coordinates in the QS simulation are only 
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FIG. 2. (color online) Rheology and velocity fluctuations at <j> = 0.836. (a) Relation between stress a and strain rate 7 obtained 
from MD simulations for £ = 10 _1 , at v = 2/3 (■), v = 1 (•) and v = 4/3 (a), and for £ = 10 -3 at 1/ = 1 (o). A small strainrate 
regime can be identified where cr = rff . The solid lines correspond to this expression, with <t] independently determined from 
the QS simulation using Eq. [2] (b-c) Probability distribution function P(A) of the rescaled velocity fluctuations A = Sv/'y at 
v = 1, for (b) C = 10 _1 and (c) \ = 10" 3 . 



available at discrete steps, one has to define an effective 
particle velocity v qs = jSf/ 70 from the particle displace- 
ment Sf during such a single step. Therefore, A is also 
the displacement rescaled by the strain interval 7 and 
P(A) is the van Hove function. Note that the viscosity 
is related to the {v + l)st moment of the velocity fluctu- 
ations, and thus of P(A) (Eq.[2]). 

Results - To characterize statistically the trajectories 
we consider the probability distribution for particle veloc- 
ities, P(A). We concentrate on the velocity component 
in the gradient direction (y-componcnt), which automati- 
cally eliminates trivial particle motion due to the average 
flow field. When the strain rate is small enough, P(A) 
reaches a limiting form (dotted line in Fig.[2)>c), which is 
directly related to the small strain-rate power-law regime 
of Fig. [2^,. Whenever the rheology a(j) deviates from this 
asymptotic behaviour, the distribution function P(A) de- 
viates from its asymptotic form as well. Interestingly, the 
approach towards this asymptotic form is rather differ- 
ent in the weakly damped, i.e. shcar-thickening, system 
(FiglJjD and • in Figl2^,) as compared to the strongly 
damped, shear-thinning system (Fig. [2k and o in Fig. [2^,). 

The velocity fluctuation PDFs obtained in the QS sim- 
ulation (solid lines in Fig [3Jd-c) arc similar to those ob- 
tained in MD in the limit of vanishing 7. In particular, 
the sharp shoulder at A w 4 is well reproduced for both 
strongly (£ = 10 _1 ) and weakly (£ = 10~ 3 ) damped sys- 
tems. Furthermore, as shown in Fig. [3l the small strain- 
rate form of P(A) only weakly depends on the value of 
the exponent v. Again, the most pronounced feature is 
the shoulder, which is nearly identical in all four simu- 
lations. However, small differences between MD and QS 
remain especially for small damping (Fig. [2t) or v > i 
(Fig. [3]). Here, a rise at small A — > is observed at the 
smallest strainrates, which is not present in the QS sim- 
ulation. Importantly, such small velocities do not con- 



tribute to the second moment of the distribution and 
therefore are irrelevant for the dissipated energy and the 
viscosity (Eq. [2]). By way of contrast, these small differ- 
ences may be important for the number of inter-particle 
contacts. As it turns out, the coordination number Z, 
which is the number of contacts per particle (taken from 
simulation snapshots), strongly varies with either v or 
C and is also different in the QS simulation; its value 
changes from Z- lso is 4, i.e. close to the isostatic state, 
down to small values Z < 1. 

We conclude that the overall features of the particle 
trajectories in the MD simulations are statistically com- 
parable to those in the QS simulation. The small strain- 
rate power-law fluid regime (Newtonian for v = 1) should 
therefore be considered as a true quasi-static limit, which 
is by no means obvious. In fact, the QS limit seems 
much better defined here ((f> < <j> c ) than in the plastic 
flow regime (0 > C ), where QS simulations have usually 
been applied, but where they suffer from a dependence 
on system size §1 [l^ ■ 

One important consequence of the equivalence MD-QS 
is that one set of QS trajectories can be used to deter- 
mine the flow rheology for different values of v. Fig. [2^ 
compares the rheology obtained using MD (data points) 
and QS simulations (solid line). They nicely collapse on 
each other when MD simulations are considered in the 
limit of small shear rate. FigfT] shows the viscosity 77 de- 
termined from both simulations, as a function of volume 
fraction (p. Beyond noting the quality of the collapse, one 
observes that the viscosity diverges with <p c — </>, with a 
scaling exponent ~ 2.2 consistent with the values mea- 
sured experimentally. 

Eq.Q, giving the power dissipated per unit volume, 
leads to the scaling law r\ ~ (<5w 1+l/ ), which connects 
the divergence of the macroscopic viscosity to the scal- 
ing law followed by the microscopic particle motion Sv ~ 
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FIG. 3. Comparison of the probability distribution function 
P(A) of the rescaled velocity fluctuations A = (v — «flow)/7 
obtained for the different computational models at different 
volume fractions cj>. Measurements are performed in the low 
shear rate asymptotic regime, for 7 = 10 -6 and £ = 0.1. The 
three values of v correspond to the dissipation mechanism of 
Eq. [1] the LF label refers to the lubrication like mechanism 
of Eq. |3] and QS, to quasi-static simulations. 



(4> c — <f>)~P ■ Thus, the seemingly harmless power balance 
turns into a relation between the exponents controlling 
the divergence of velocity fluctuations and that of vis- 
cosity: a = (3(1 + v). We have recently shown that 
(3 ~ 1.1 [2^], which gives (for v = 1) a c± 2.2, con- 
sistent with the exponent extracted from the MD data 
in Fig. [T] Note however, that subdominant corrections 
can lead to apparent exponents that, in the considered 
range of densities, may not reflect the true asymptotic 
behavior I . 

Discussion - The small strain-rate rheology, cr = rff v , 
as well as the divergence of the viscosity, r\ ~ 8(f/^ 1+v ' , 
depend on the value of v. On the other hand, the under- 
lying particle trajectories are hardly affected by changing 
v. This points to a certain decoupling between parti- 
cle trajectories and dissipative process. In this picture, 
the statistical properties of trajectories are largely gov- 
erned by the structural singularity of random close pack- 
ing and the lack of space available for particle motion. 
On the other hand, system-specific dissipation mecha- 
nisms affect the rheological properties via the dissipated 
energy along these geometrically predetermined trajec- 
tories. Certainly, such a decoupling cannot be realized 
in a perfect manner, as shown by the small differences 
of P(A) at small A and in the tails (Fig. [2] and Fig. EJ. 
Nevertheless, it seems to be strong enough such that the 
various flow curves (Fig. [2k) and the viscosity (Fig. [J) 
can accurately be predicted from the sole knowledge of 
one set of QS trajectories. The scaling law relating the 
viscosity to the volume fraction is a directly testable pre- 



diction of the central idea of this letter. It suggests to 
measure the rheology of particles suspended in a non- 
Newtonian solvent like a polymer melt or a visco-plastic 
fluid. 

In order to investigate the universality of the decou- 
pling phenomenon, we have conducted additional simu- 
lations with a damping that describes a modified lubri- 
cation force 27J between neighboring particles i and j: 



F diss (^) = -Cn 4 



{Vi -Vj)] 



(3) 



In agreement with previous (overdamped) simulations [H, 
[iH, the spatial velocity-correlation function is qualita- 
tively different in the models of Eqs.([T]) and ((3]) (not 
shown). Still, the probability distribution P(A) (FigJ3]) is 
remarkably similar across all models considered. In par- 
ticular, the scale of velocity fluctuations increases on ap- 
proaching the critical volume-fraction <p c , and the overall 
agreement between the different curves seems to improve. 
This supports our interpretation of the role of close pack- 
ing for the particle trajectories. Furthermore, starting 
with the QS trajectories, a calculation similar to Eq. ([2]) 
can again be used to predict the viscosity. As FigQ] il- 
lustrates (open circles) this calculation is quite accurate 
but slightly overestimates the true viscosity (triangles), 
roughly by a factor wl.5. 

In conclusion, singular velocity fluctuations cause a dy- 
namical contribution to the di verg ence of the viscosity, 
as independently noted in Ref. [18(. These velocity fluc- 
tuations are surprisingly conserved across different com- 
putational models, which we explain with geometric fea- 
tures and the lack of available space close to the jamming 
transition. Our results complement those obtained in 



Ref. |18[, for frictionlcss hard spheres. In that system, the 
jamming transition coincides with the isostatic threshold 
and the flow properties can be related to the geometry 
of the contact network. The coordination number Z is 
then the relevant control parameter. Our results bridge 
the gap between such an ideal system and those where 
inertia and elasticity lead, for a given volume fraction </>, 
to strong changes of Z . Indeed, in our simulations, the 
value of Z has no immediate predictability and <f> is the 
only relevant parameter. Our results open the promising 
perspective of the existence of an inherent contact net- 
work which would govern the topography of the energy 
landscape. Such a concept would establish the missing 
connection between volume fraction and connectivity. 
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